{
    "cells": [
        {
            "cell_type": "markdown",
            "id": "5be6c57b-7cae-4fc1-b78f-899becabc6ee",
            "metadata": {},
            "source": [
                "# TorchSWE case study\n",
                "\n",
                "\n",
                "[TorchSWE](https://github.com/piyueh/TorchSWE) is a shallow-water solver created by Dr. Pi-Yueh Chuang and Prof. Lorena Barba that solves the vertically averaged Navier-Stokes equations using MPI and CuPy. It can simulate free-surface water flow in rivers, channels, and coastal areas, as well as model flood inundation. Given a topography, TorchSWE can predict flood-prone areas and the height of water inundation, making it a valuable tool for risk mapping.\n",
                "\n",
                "High-resolution numerical simulations\u2014such as those on real topographies requiring hundreds of millions of data points\u2014demand distributed computation across multiple GPUs. Although scalability is achievable with MPI4Py and CuPy, this approach requires manually partitioning the problem and managing inter-GPU data communication, which are complex and error-prone tasks.\n",
                "\n",
                "cuPyNumeric enables a distributed implementation of TorchSWE using only NumPy operations, without the complexities of MPI+CuPy. After porting TorchSWE to cuPyNumeric by removing all domain decomposition logic, it scaled effortlessly across multiple GPUs and nodes without further code modifications. This scalability enabled high-fidelity simulations exceeding 1.2 billion data points using 32 GPUs, allowing researchers to tackle critical scientific problems in flood inundation modeling without needing specialized distributed computing expertise. Overall, the cuPyNumeric implementation reduced the lines of code by over 20%, and simplified development and maintenance by eliminating complex logic for managing distribution and communication.\n"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "0402fb01-748b-48d9-9caa-80e7510ade80",
            "metadata": {},
            "source": [
                "\n",
                "<h2>Deep dive into the TorchSWE code implementation</h2>\n",
                "\n",
                "<h3> Original code details</h3>\n",
                "\n",
                "TorchSWE uses stencil operations to model shallow-water equations on a 2D grid, where each point is updated based on neighboring values, simulating water flow dynamics. The stencil computations are structured to update each grid cell iteratively, based on data from surrounding cells, mimicking fluid behavior over time. Below is an example that mimics the basic structure of the stencil logic from the TorchSWE repository:\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 1,
            "id": "640f0b62-f70f-4d8a-86c5-7b4739e60a33",
            "metadata": {},
            "outputs": [],
            "source": [
                "import numpy as np\n",
                "\n",
                "# Example dimensions for the grid\n",
                "nx, ny = 128, 128\n",
                "grid = np.ones((nx, ny))  # Initialize the grid with \"1\"\n",
                "\n",
                "# Stencil operation\n",
                "for i in range(1, nx - 1):\n",
                "    for j in range(1, ny - 1):\n",
                "        grid[i, j] = (\n",
                "            grid[i + 1, j] + grid[i - 1, j] + grid[i, j + 1] + grid[i, j - 1]\n",
                "        ) / 4"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "0281b3f4-5a48-40cc-9ec8-0fc9d7fd760c",
            "metadata": {},
            "source": [
                "This code iteratively updates cell `h[i, j]` using adjacent cells, representing a basic averaging stencil operation that can be extended to various boundary conditions and flow dynamics in the shallow-water model. For full context, refer to [TorchSWE on GitHub](https://github.com/piyueh/TorchSWE).\n",
                "\n",
                "Parallelizing stencil operations for multi-GPU systems is challenging. When arrays are partitioned across multiple GPUs, any update to a cell requires the updated values to be shared between GPUs to maintain consistency across boundaries. This communication overhead and synchronization make parallelizing stencil code complex and difficult to implement efficiently on multi-GPU architectures.\n",
                "\n",
                "Below, we outline TorchSWE\u2019s MPI4Py logic in more detail  to highlight the complexity involved in this implementation.\n",
                "Here\u2019s an example code snippet that mirrors the TorchSWE MPI logic, implementing a simple MPI stencil operation from above:\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 2,
            "id": "0d7db631-3ae9-41ca-a0f1-07390349fbd0",
            "metadata": {},
            "outputs": [],
            "source": [
                "from mpi4py import MPI\n",
                "import cupy as cp\n",
                "\n",
                "num_timesteps = 10\n",
                "\n",
                "\n",
                "def set_device(comm: MPI.Comm):\n",
                "    # Device selection for each rank on multi-GPU nodes (TorchSWE-specific)\n",
                "    n_gpus = cp.cuda.runtime.getDeviceCount()\n",
                "    local_rank = comm.Get_rank() % n_gpus\n",
                "    cp.cuda.runtime.setDevice(local_rank)\n",
                "\n",
                "\n",
                "comm = MPI.COMM_WORLD\n",
                "rank = comm.Get_rank()\n",
                "size = comm.Get_size()\n",
                "\n",
                "# Determine grid size and decompose domain\n",
                "gnx, gny = 126, 126  # global grid dimensions\n",
                "local_nx, local_ny = gnx // size, gny  # local grid dimensions per rank\n",
                "local_grid = cp.ones((local_nx + 2, local_ny + 2))  # with halo boundaries\n",
                "\n",
                "# Set up MPI data types and boundaries\n",
                "send_type, recv_type = (\n",
                "    MPI.DOUBLE.Create_subarray(\n",
                "        (local_nx + 2, local_ny + 2), (local_nx, local_ny), (1, 1)\n",
                "    ),\n",
                "    MPI.DOUBLE.Create_subarray(\n",
                "        (local_nx + 2, local_ny + 2), (local_nx, local_ny), (1, 1)\n",
                "    ),\n",
                ")\n",
                "send_type.Commit()\n",
                "recv_type.Commit()\n",
                "\n",
                "# Stencil computation loop\n",
                "for timestep in range(num_timesteps):\n",
                "    # Boundary exchange with non-blocking sends/receives\n",
                "    reqs = []\n",
                "    if rank > 0:\n",
                "        reqs.append(comm.Isend(local_grid[1, :], dest=rank - 1))\n",
                "        reqs.append(comm.Irecv(local_grid[0, :], source=rank - 1))\n",
                "    if rank < size - 1:\n",
                "        reqs.append(comm.Isend(local_grid[local_nx, :], dest=rank + 1))\n",
                "        reqs.append(comm.Irecv(local_grid[local_nx + 1, :], source=rank + 1))\n",
                "\n",
                "    # Ensure all sends/receives are complete\n",
                "    MPI.Request.Waitall(reqs)\n",
                "\n",
                "    # Perform stencil operation\n",
                "    for i in range(1, local_nx + 1):\n",
                "        for j in range(1, local_ny + 1):\n",
                "            local_grid[i, j] = 0.25 * (\n",
                "                local_grid[i - 1, j]\n",
                "                + local_grid[i + 1, j]\n",
                "                + local_grid[i, j - 1]\n",
                "                + local_grid[i, j + 1]\n",
                "            )\n",
                "\n",
                "# Clean up MPI data types\n",
                "send_type.Free()\n",
                "recv_type.Free()\n",
                "MPI.Finalize()"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "660621f9-2bc9-49a3-be59-cde1ce87df65",
            "metadata": {},
            "source": [
                "This example follows TorchSWE's approach to domain decomposition and parallelization as in the original implementation. It starts with MPI initialization and sets up logic to manage GPU assignment per rank, dividing the global grid into subdomains. Each rank is responsible for a local subgrid with added halo rows to hold neighboring data. Once the domain is decomposed, the user must ensure proper communication of data at processor boundaries, accounting for datatype differences between CuPy and MPI4Py. For optimal performance, the appropriate type of point-to-point communication, such as non-blocking send/recv, must be selected, as incorrect implementation can cause deadlock. Users must also handle varying numbers of neighboring ranks on domain boundaries and ensure data exchange across mesh, topography, and solution variables. Non-blocking `Isend` and `Irecv` functions handle boundary data exchanges, allowing each rank to receive necessary data for stencil computations. After a `Waitall` synchronization step, each rank performs computations on its subdomain. Finally, custom MPI data types are freed, and `MPI_Finalize()` concludes the environment.\n",
                "\n",
                "The actual TorchSWE code has additional complexities specific to its use of multiple arrays, GPU memory management, one-sided communications etc.\n",
                "For the complete implementation, you can refer to the [TorchSWE repository](https://github.com/piyueh/TorchSWE).\n",
                "\n",
                "Explicit distributed logic, like that in TorchSWE, is difficult to debug and maintain throughout the lifespan of simulation codes. Most applications, including TorchSWE,  require specialized validation tests to ensure correct outputs. This results in significant programming effort and further complicates development. \n"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "e93aa24e-fc18-4f69-819d-59b5997aa087",
            "metadata": {},
            "source": [
                "<h3>cuPyNumeric Implementation</h3>\n",
                "\n",
                "In the [cuPyNumeric version of TorchSWE](https://github.com/shriram-jagan/TorchSWE), stencil operations are implemented using distributed array handling from cuPyNumeric, simplifying the code and removing the need for manual partitioning or boundary synchronization. The code operates similarly to NumPy slicing but scales across multiple GPUs. For example, the stencil computation in this version would typically involve using simple array slices like below (instead of the nested loops with integrated MPI logic as in the original implementation).\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "id": "b6e15757-a681-4a09-9f82-6304adf82fb4",
            "metadata": {},
            "outputs": [],
            "source": [
                "import cupynumeric as np\n",
                "\n",
                "# Example dimensions\n",
                "nx, ny = 128, 128\n",
                "\n",
                "# Initialize the array h\n",
                "grid = np.ones((nx, ny))\n",
                "\n",
                "# Stencil operation using slicing\n",
                "grid[1:-1, 1:-1] = (\n",
                "    (\n",
                "        grid[2:, 1:-1]  # Below\n",
                "        + grid[:-2, 1:-1]  # Above\n",
                "        + grid[1:-1, 2:]  # Right\n",
                "        + grid[1:-1, :-2]  # Left\n",
                "    )\n",
                "    / 4\n",
                ")"
            ]
        },
        {
            "cell_type": "markdown",
            "id": "f29f5387-3408-4bff-948d-55519412de31",
            "metadata": {},
            "source": [
                "This operation is automatically managed across nodes and GPUs without needing MPI-specific code. More details can be found in the [cuPyNumeric port of TorchSWE](https://github.com/shriram-jagan/TorchSWE).\n",
                "\n",
                "The cuPyNumeric version of TorchSWE eliminates 600 lines of code related to domain decomposition, communication, synchronization, and validation that would otherwise be needed when using MPI4Py with CuPy. These 600 lines require substantial knowledge of distributed computing from domain scientists. By using cuPyNumeric, the simplified NumPy code scales efficiently to 1024 GPUs, making high-fidelity flood modeling accessible without requiring specialized expertise in distributed systems."
            ]
        },
        {
            "cell_type": "markdown",
            "id": "7e5d6565-ceda-4b61-8826-b6ae5aff3c83",
            "metadata": {},
            "source": [
                "<h2>Conclusion</h2>\n",
                "\n",
                "cuPyNumeric significantly simplifies the development and maintenance of distributed simulations, such as TorchSWE, by abstracting complex parallelization, synchronization, and communication logic. This eliminates the need for specialized HPC knowledge and reduces the risk of errors, allowing domain scientists to focus on their research. With cuPyNumeric, large-scale simulations can scale efficiently across large HPC systems, enhancing productivity, reducing programming effort, and lowering development costs. \n",
                "\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "id": "eb3a186a-3ea7-4150-8ec0-7760ad2adf1f",
            "metadata": {},
            "outputs": [],
            "source": []
        }
    ],
    "metadata": {
        "kernelspec": {
            "display_name": "Python 3 (ipykernel)",
            "language": "python",
            "name": "python3"
        },
        "language_info": {
            "codemirror_mode": {
                "name": "ipython",
                "version": 3
            },
            "file_extension": ".py",
            "mimetype": "text/x-python",
            "name": "python",
            "nbconvert_exporter": "python",
            "pygments_lexer": "ipython3",
            "version": "3.12.7"
        }
    },
    "nbformat": 4,
    "nbformat_minor": 5
}
